@include "globals.awk"
@include "lib.awk"
@include "readcsv.awk"
@include "dist.awk"
@include "tile.awk"

function _cube(  o,what,who,_Cube) {
  what="CUBE: map data to N-dimensional hypercube"
  who ="(c) 2012 tim@menzies.us"
  if (args("--c,0,--help,0,-d,data/nasa93.csv,"\
           "--dump,0,-seed,1,-p,10,--n,1,-m,3,-r,-1,-k,1",o)) {
    if(o["--help"]) return prints(what,who," ",
     " --c       copyright information",
     " -d FILE   data file",
     " --help    print help",
     " --n       disable normalization",
     " -m N      max dimensions; default=3",
     " -p N      probes for dimension#1; default=1",
     " -seed N   random number seed; default=1")
    if (o["--c"])
      return copyleft(what,who)
    seed(o[ "-seed"])
    probes   = o["-p"]
    normp    = o["--n"]
    maxSides = o["-m"] + 1
    rnn      = o["-r"]
    dump     = o["--dump"]
    k        = o["-k"]
    cube(o["-d"],_Cube)
}}
function cube(f,   _Cube,_Rows,    i,rows,cols,hi1,lo1,out,chops) {
  makeCube(f,_Cube,_Rows)
  lessRows(_Cube,_Rows,rows)
  print " :r ",rnn," :l ",length(rows)
  lessCols(_Cube,_Rows,cols,hi1,lo1)
  for(i in all)
    out[i]=cubeknn(all[i],_Cube,_Rows,rows,cols,hi1,lo1) 
  pairs("0.1,-,0.3, ,0.5, ,0.7,-,0.9, ",chops)
  print tile(out,chops,30,"%3s",0,1200)
}
function cubeknn(row1,_Cube,_Rows,rows,cols,hi1,lo1,\
                 klass,kth,i,got,want,g) {
  g = klass1(gp)
  cubekth(row1,_Cube,_Rows,rows,cols,hi1,lo1,kth)
  for(i in kth)
    got += d[kth[i]][g]
  got  = got/k
  want = d[row1][g]
  return int(100*abs(want - got)/want)
}
function cubekth(row1,_Cube,_Rows,rows,cols,hi1,lo1,kth,\
                 i,row2,n,diff,sort) {
  for(i in rows) {
      row2 = rows[i]
      if(row1==row2) continue
      diff[i]["="] = row2
      diff[i]["x"] = euclidean(cols[row1],cols[row2],hi1,lo1)
  }
  n = asort(diff,sort,"xsort")
  k = n < k ? n : k
  for(i=1;i<=k;i++)
    kth[i] = sort[i]["="]
}
function euclidean(a1,a2,hi1,lo1,    i,d,one,two,n) {
 for(i in a1) {
   n++
   one = a1[i]; one = (one - lo1[i])/(hi1[i] - lo1[i] + 0.0001)
   two = a2[i]; two = (two - lo1[i])/(hi1[i] - lo1[i] + 0.0001)
    d += (one - two)^2
  }
  return sqrt(d)/sqrt(n)
}
function makeCube(f, _Cube,_Rows,    p) {
  readcsv("cat " f,_Rows)
  p         = dim0(_Cube,_Rows)
  corners[1] = lhs[p]
  corners[2] = rhs[p]
  anotherSide(probes, rhs[p],_Cube,_Rows)
}
function lessRows(_Cube,_Rows,selected,\
                           p,i,tmp) {
  for(p in x) 
     rnns(sorted[p],tmp)
  for(i in tmp)
    if (tmp[i] > rnn)
      selected[i]= i
}
function  rnns(x,count,   \
              i,j,k,down,here,up,ix,jx,kx,left,right,max) {
  max = length(x)
  for(j=2;j<max;j++)
    rnn(count,
        x[j-1]["x"],x[j]["x"],x[j+1]["x"],
        x[j-1]["="],x[j]["="],x[j+1]["="])
}
function rnn(count,down,here,up,ix,jx,kx,\
             left,right) {
    left  = here - down
    right = up - here
    if (left > right) 
      count[kx]++ 
    else if (left < right) 
      count[ix]++ 
    else { 
      count[ix]++
      count[kx]++ }
}
function lessCols(_Cube,_Rows,projected,hi1,lo1,  \
                  g,i,j,one,p,q,line,first,n,order,val) {
  n=sides(x,order)
  g = klass1(gp)
  for(i in all) {
    one=all[i]
    if (dump) line=""
    if(!first) {
      first = "#"
      for(q=1;q<=n;q++) {
        p = order[q]
        hi1[p] = -1 * ( lo1[p] = 10**32 )
        first = first "$side_"p OFS }
      if (dump)
        print first "$!class"
    }
    for(q=1;q<=n;q++) {
      p = order[q]
      val = x[p][one]
      projected[one][p] = val
      if (val < lo1[p]) lo1[p]=val
      if (val > hi1[p]) hi1[p] = val
      if (dump)
        line = line sprintf("%6.4f",val) OFS
    }
    if (dump)
      print line d[one][g] "# " one | "sort -n -k 2"
  }
}
function sides(x,order,    tmp,i) {
  for(i in x) tmp[i]=i
  return asort(tmp,order)
}
function anotherSide(p,old,_Cube,_Rows,  new) {
  if (length(corners) >= maxSides) 
    return p 
  new = furthestFromAllCorners(_Cube,_Rows)
  if (!new) return p
  push(new,corners)
  p++
  pole(p,old,new, _Cube,_Rows)
  return anotherSide(p,new,_Cube,_Rows)
}
function furthestFromAllCorners(_Cube,_Rows,   \
                                i,one,now,most,out) {
  for(i in all) {
    one = all[i]
    if (has(corners,one))
      continue
    now = fromAllCorners(one,_Cube,_Rows)
    if (now > most) { out = one; most = now}
  }
  return out
}
function fromAllCorners(here,_Cube,_Rows,   there,out) {
  for(there in corners)
    out += Dist(here,there, _Cube,_Rows)
  return out
}
function dim0(_Cube,_Rows,\
              best,p,this,out) {
  best = 10**32
  for(p=1;p<=probes;p++) {
    this = anyPole(p,_Cube,_Rows)
    if (this < best) { 
      best= this; 
      out= p 
    } else 
      oneLess_Cube(p,_Cube)
  }
  for(p=1;p<=probes;p++) 
    if (err[p] != best)
      oneLess_Cube(p,_Cube)
  return out
}
function show_Pole(p,_Cube) {
    print "\n"p
    o(x[p],"x")
    print ":lhs " lhs[p] " :rhs " rhs[p]
}
function oneLess_Cube(p,_Cube) {
  delete poles[p]; delete err[p]
  delete     a[p]; delete   b[p]; delete c[p]
  delete     x[p]; delete   y[p]
  delete     l[p]; delete   r[p]
  delete   lhs[p]; delete rhs[p]
  delete sorted[p]; delete corners[p]
}
function anyPole(p,_Cube,_Rows,   one,two) {
     one  = any(all)
  do {two = any(all)} while (one == two)
  return pole(p,one,two,_Cube,_Rows)
}
function pole(p,one,two,_Cube,_Rows,\
              most,least,i,three,aa,bb,cc,xx,yy) {
  l[p] = one
  r[p] = two
  c[p] = cc = Dist(l[p],r[p],_Cube,_Rows)
  most = -1 * (least = 10**32)
  for(i in all) {
    three = all[i]
    if ( three == l[p] ) {
      aa = xx = 0; bb = cc
    } else {
      if (three == r[p] ) {
        aa = xx = cc; bb = 0
      } else {
        aa = Dist(three, l[p], _Cube, _Rows)
        bb = Dist(three, r[p], _Cube, _Rows)
        xx = (aa^2 + cc^2 - bb^2)/(2*cc + 0.0000001)
      }}
    yy = sqrt(abs(aa^2 - xx^2))    
    x[p][i]  = xx
    sorted[p][i]["="] = three
    sorted[p][i]["x"] = xx
    if (xx < least) { lhs[p] = three; least = xx} 
    if (xx > most ) { rhs[p] = three; most  = xx}
    y[p][i]  = yy
    err[p]  += yy/length(all)
  }
  asort(sorted[p],sorted[p],"xsort")
  return err[p]
}
function Dist(i,j,_Cube,_Rows,   out) {
  if ( (i,j) in cache ) 
    out = cache[i,j]
  else {
    out = dist(i,j,_Rows,normp)
    cache[i,j] = cache[j,i] = out
  }
  return out
}
